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Abstract 

We present results of simulations of the quenched quarkonium spectrum at 
two values of the lattice spacing and for quark masses around rUc using the 
tadpole improved clover action. Attention is focussed on the lowest lying S 
and P states, and the triplet fine structure is obtained for the first time using 
a relativistic action. 
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I. INTRODUCTION 



Lattice QCD simulations provide a systematically improvable method for predicting low 
energy phenomenology directly from QCD, providing a means to verify that QCD correctly 
reproduces experiment. We are primarily interested here in the simulation of charmonium, 
which remains a challenge for lattice QCD since the quark mass is sufficiently large to be 
difficult to resolve on coarse lattices, while the system is rather relativistic, causing problems 
for the NRQCD approach. 

The angular momentum states calculated here are the ^5*0, ^Si, ^Pq, ^Pi, and ^P2 of 
lowest radial excitation. This is more extensive than typical lattice calculations using a 
relativistic action, which have generally been limited to the ^5*0, ^Si and spin-averaged S-P 
splittings though Bhattacharya et al produce a value for the Xci — Xco fine splitting 

using unimproved Wilson quarks. 



A. Simulation Techniques for Charm and Beauty 

A discrete action approximating the continuum action is used, where finite differences 
are substituted for derivatives, summations for integrals and so on, giving the 0{a) improved 
Wilson quark action 

Si.u = ^T.^i^)M{x,ymy), (1) 



M{x,y) 



1 - 



a F 



(2) 



where k = 2moa+8 ' g^iuge group elements U^{x) = e^<^9Ati.(x+^fi) ^ expansion of the 

discrete action in terms of continuum fields and their derivatives gives a natural classification 
of errors in powers of the lattice spacing a 



latt 



5, 



QCD 



ij{x),ij{x),A^{x) +aSi tP{x),ij{x),A^{x) +0{a^). (3) 



We have chosen a lattice action which reduces to QCD as a — 0, and may eliminate the 
leading errors due to by adding appropriate (higher dimension) operators. The clover 
term 



vCsw—(T^ivF^,ui). 



(4) 



gives an 0(a) contribution to the action so is irrelevant to the continuum limit; its coefficient 
may be selected to minimise the leading discretisation effects contained in Si. The fixing of 
the clover coefficient has been performed in various ways - to different orders in perturbation 
theory using mean- field improvement 0, and most recently using non-perturbative 

improvement criteria |T0[ , giving a much improved continuum scaling behaviour of the light 
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hadron sector |TT]]. This categorisation of errors in powers of the lattice spacing a, however, 
introduces an upper hmit for quark masses that can be accurately simulated, amq <^ 0(1). 
Within the quenched approximation contemporary simulations typically have amcharm > |; 
leaving them apparently susceptible to significant discretisation effects, which can only be 
removed by brute force extrapolation to the continuum limit. In recent years the problem 
of simulating heavy quarks in lattice QCD has therefore been tackled by discretising the 
effective theory, NRQCD - essentially performing a double (power counting) expansion in 
both the lattice spacing and in \ |jl2|, producing very wide ranging spectrum predictions for 
quarkonia [|T^-|T5[| . However, the non-relativistic approach has recently been demonstrated 
to suffer from significant relativistic corrections for charmonium , leading to inconsistent 
results when comparing different orders of the nonrelativistic action fl^, while the relativistic 



corrections to T appear to be well under control. 

If we make a power counting expansion, in the manner of Fermilab of the improved 
Wilson action in both and a simultaneously, however, we can see that the situation 
for heavy valence quarks is not as bad as naively supposed. Consider the tree-level low 
momentum expansion of the energy, (while this applies to the free case, the errors at large 
rrioa are certainly indicative of errors in the interacting theory). 



where 

ami = log (1 + amo), (6) 

am, = e-^ sinh(amO 
1 + sinh(ami) 

„^ sinh(ami) , , 

arriB = e""^' — ^— ^ r (8) 

^ + (-'SW sinh(ami) 

am, = e-^ '^^^ , (9) 



1 + e"™i sinh(ami)(l + sinh (ami)) + 26^"™! sinh(ami) 
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where the A^^ term violates rotational symmetry at O(^) and B^^ contains order 
spin-dependent corrections. 

In this light, the accuracy of the action depends on both the size of am and on the 
size of f ^. In a sufficiently non-relativistic system the above Hamiltonian will, at tree level, 
accurately describe a particle of mass m2 = ms- If am <^ 1 the action is approximately 
Lorentz symmetric [] and the coefficients of higher order terms in the above non-relativistic 



^assuming that aAqcD 1 
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action will accurately match those of the continuum expansion. The essential point for 
charm however is that we have to estimate the size of errors in the subleading terms relative 
the terms containing m2 and m^. The errors will depend on both the difference between 
the higher order coefficient from that of the continuum Hamiltonian for a particle of mass 
777-2, and on the size of an 0(f^") prefactor. At tree level, therefore, we can estimate the 
fractional error introduced in the kinetic energy due to discretisation effects in the action as 
follows. 

Akin - — — [i-yj) 
Taking a charm quark mass of around 1.5 GeV and an estimate for (f^) ~ 0.3 ""^^l^™ 



(which approximates the data in ||T8[ with sufficient accuracy to satisfy our needs), we plot 
the fractional error in the kinetic energy. Akin as a function of am for both of the lattice 
spacings we use in figure |I[ The leading order error in the kinetic energy is estimated to 
be between 5% — 10%. We note that had we used the pole mass, the leading correction 
would have arisen from the mi — m2 discrepancy rather than the m2 — m^ discrepancy, and 
therefore would not have been suppressed by f^. Both analyses would, of course, have given 
the same continuum limit; however the use of 7772 will for this reason have a milder lattice 
spacing dependence. 

The disadvantage of any simulation involving Wilson quarks with amo > 0(1) is that 
the a dependence of observables is not likely to have a simple low order polynomial form. 
There is no hierarchy in the orders of a due to terms (ma)"', which is manifested in 
the non-linear dependence of the various mass definitions on the bare quark mass. These 
different mass definitions rapidly become identical as ma is reduced below 1, so that, for 
charm, studying the a dependence is sensible with current lattice spacings. For bottomonium 
results, however, amQ > 1 in both simulations, and we must rely on the power counting 
estimates of the correction terms to quantify errors. Continuum scaling plots provide only 
circumstantial evidence that the systematic errors in the bottom simulations are controlled. 
A more rigorous approach would be to form the union of bounds estimated for different 
error terms on each of the lattices; in this respect simulations in which amQ > 1 must be 
treated in a similar way to effective theories - the fact that the continuum limit exists for 
the Fermilab approach [|l| is not in practice an advantage for bottom simulations since we 
do not know how to extrapolate /ro?77 this regime. 

Ideally we would make a corresponding estimate for the errors in the spin-dependent 
pieces of the action, but the subleading spin- dependent term for the Wilson action has not 
presently been calculated. However, it is thought that a similar size systematic error is likely 
in the spin-dependent Hamiltonian. 

Some progress has recently been made in treating Wilson type quarks to one loop in 
perturbation theory for arbitrary quark mass. In particular Kronfeld et al [|l9l have calcu- 



lated the 1-loop renormalisation of the terms 7r7i and 7772, showing that the discrepancy is 



^we are dealing with on-shell correlation functions of valence quarks 
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dominated by the tree level result. The Kuramashi pO| calculation gave bilinear current 
renormalisation constants for Wilson quarks, and somewhat interestingly, noted that it was 
necessary to match m2 to the continuum pole mass in order to obtain the same IR behaviour 
for the vertex and wavefunction in the lattice theory as in the continuum theory. In order 
to obtain more rigorous control of the errors in charmonium simulations it would, of course, 
be interesting to have a one-loop calcution for m^, and (at least) a tree level estimate of the 
subleading spin-dependent pieces for Wilson type actions. 

In the remainder of this paper we present a study of quarkonium systems using the 
tadpole improved clover action |^,[^ at both /3 = 6.0 and [3 = 6.2, corresponding to a ~ 0.1 
fm and a ^ 0.07 fm, analysed in the manner suggested by ||T|. In the calculation at /3 = 6.0 
similar operators to those used in NRQCD calculations are applied enabling us to produce a 
wider ranging spectrum than previous relativistic calculations. The scaling behaviour of the 
hyperfine splitting between the two lattices is discussed, and estimates of the consistency 
of continuum quenched results with experiment are made. We note that while El-Khadra, 
Kronfeld, Mackenzie and others ||l],|19[ propose improving the action for on shell quantities 



by breaking 0(4) rotational symmetry in the Fermilab action, their simulations to date have 
involved the same re-interpretation of the standard relativistic quark action as we use here, 
and the results are therefore directly comparable. 

II. QUENCHED QUARKONIUM SPECTRUM AT (3 = 6.0 
A. Simulation Details - f3 = 6.0 calculation 

The first simulation was performed using 499 sample gauge configurations from a 
quenched distribution at (3 = 6.0 on a 16^ x 48 lattice. Five heavy quark masses were 
used, as given in Table |I|, with the tadpole improved clover action with Csw = 1.47852, 
corresponding to uq = 0.8778 taken from the plaquette. 

The timesliced residue 

E I E.',t' M{x, t; x', t')^{x', t') - 7]{x, t) |2 

I |2 _ x'gL^ (11') 

where the normalisation of the metric reflects the exponential fall off of the pseudoscalar 
correlator, was used to define the convergence criterion. Using the variance of previous sim- 
ulation data to estimate the expected statistical error in the correlation functions, we found 
that requiring the timesliced residue on timeslice 24 ( the noisiest ) of rt=24 < O(10~^^), gives 
an error due to convergence when comparing to an extremely tightly converged propagator 
that is orders of magnitude below the expected variance. This was subsequently used as a 
convergence criterion and corresponds to a traditional residue of between 10~® and 10~^^ 
depending on the quark mass, significantly tighter than the 10"^ obtainable using 32 bit 
precision. 

Both local and fuzzed sources were used . Fuzzed smearing uses Michael- Teper fuzzed 



links to covariantly transport the quark fields n links along each of the axes to form a face- 
centred cubic smearing function. The fuzzing radius n = 4 was selected to optimise the 
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plateau for P-states using a small number of configurations, which, in combination with 
sink fuzzing, allowed a 2 x 2 smearing matrix to be constructed. 



in n\-( {OL{x)OLm {OAx)OLm \ 

Three additional propagators, from covariant derivative sources 

M^) = E [H^n^ + ly)- Ulix - mx - I y)] 6{y, 0), (13) 
y 

were generated at k = 0.12600, lying closest to charm, allowing the differential operators in 
Table || to be constructed for mesons composed of this quark and one at another value of 
K. The most important result of this is that we now have access to the angular momentum 
2 P-state via the 2~^~^ operator, in both the T and E representations. 



B. Mass Fitting of Correlators 

In what follows, all fits to correlators were made using correlated fits with the Marquardt 
Levenberg algorithm, and the bootstrap algorithm was used to estimate the errors, using 
1000 bootstrap samples. Single exponential, row fit and matrix fit models were used. 

The fit ranges were varied with tmax holding typically three values selected by inspection 
of the effective mass plots, and tmin varied over the entire range in which the Marquardt 
algorithm converged for the fit model. This procedure was carried out for each operator 
O G {0 ^, 1 , 0"*""^, l"*"*", 1"*""} and, where the data existed, for the 2"*"+ operator. 

The information was used to select an optimal fitting range by considering the stability of 
the central value and error with varying tmin, requiring Q > 0(0. 1) and x^/ dof < 0(1.5), and 
considering the consistency in the fitted masses between the double and single exponential 
fit methods. 

The simultaneous double exponential fit to the fuzzed-source local-sink, and the local- 



local correlator was the selected method, and the ranges in Table |T| were deemed optimal. 
These were subsequently used to obtain the best fit masses in Table |I^ and Table 0. The 
excited state masses given in these tables were not found to be sufficiently independent of 
the fit range to produce confident spectrum predictions. The effective mass plots obtained 
for the various angular momentum states of the ki = K2 = 0.12600 system are displayed in 
Figure |^ to Figure ^, where the open circles are the points from the local-local correlation 
function, and the fancy diamonds from the fuzzed-local. They are considered representative 
of the plots for other k, combinations. 

Since the P-states have a large splitting from the pseudoscalar, the noise grows rapidly 



with time [£2], and it is essential to perform part of the fit at early timeslices to obtain 
a good signal, requiring a well chosen fuzzing or smearing radius. Figure ^ clearly shows 
that we have resolved the complete fine structure of the x states for the first time with 
a relativistic action, since the 68% error bounds do not overlap. We, of course, obtain a 
reduced error by taking the bootstrapped difference in the fitted mass values. 
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C. Kinetic Masses 



The correlation functions were computed at various momenta allowing us to compute 
the kinetic mass by fitting the dispersion relation 

2 

Ei/) = mi + + Cp\ (14) 

The non-zero momentum behaviour of the (2600, 2600) pseudoscalar, Figure |^ is consid- 
ered typical, and the fitted lattice kinetic masses are given in Table VII| . 



D. Scale from Quarkonium S-P splitting 

When quarkonia have been simulated (using either NRQCD or the heavy Wilson ap- 
proach) quenched calculations have obtained a large slope in the -g dependence of the 
spin-averaged IS* — IP splitting while the experimental splitting in T and J/ip is thought 
to differ by only a few MeV (with reasonable assumptions for the i^h mass). The slope cal- 
culated is a quenching effect il^. It has been argued that the source of this slope is the 



incorrect running of the coupling in the quenched approximation P3| , |2^ and that the defi- 
nition of appropriate for quarkonia is that taken from the S-P splitting, in an attempt 
to eliminate quenching effects inherent in setting the scale from a low momentum quantity. 
This definition of the lattice spacing therefore depends on the mass of the heavy quarks 
used to determine it. The mass dependence of the quarkonium IP — IS* splitting and the 
corresponding inverse lattice spacing obtained in this simulation is shown in Figure |[ The 
inverse lattice spacings at /5 = 6.0 obtained for J/ip and T are 2.17(6) GeV and 2.4(1) GeV 
respectively. The value obtained using NRQCD on the same lattice was 2.5(1) for the T 
system [jl3[. The inverse lattice spacings obtained from Mp ( with the tadpole improved 



action ^ ) and ^ are 1.96(4) GeV and 2.02(3) GeV respectively. 



It has recently been suggested [113 based on potential models that the errors in the 



quenched quarkonium spectrum are dominated by an overestimation of the IS" state rel- 
ative to the physically larger states due to underestimation of the Coulomb coefficient in 
the quenched potential. As a result, even under a redefinition of the scale, the quenched 
quarkonium spectrum is inconsistent with experiment. For these reasons we analyse the 
13 = 6.0 simulation in three ways, using the p mass, the string tension, and the quarkonium 
S-P splitting to set the scale, the latter being included largely to allow comparison to other 
calculations. Finally we note that the sensitivity of both the hyperfine splitting and the 
S — P splitting to quenching makes them excellent candidates for observing unquenching 
effects in hadrons in exploratory dynamical simulations. 

E. Charm and Beauty Systems 

The extrapolations in inverse heavy quark mass for the heavy-heavy systems were found 
to be amenable to a linear fit. The hyperfine splitting, S-P splitting and 2~^^ - 1++ - 0+^ 
fine structure were calculated and extrapolated to match the kinetic spin-averaged S-state 
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mass, M2, to the physical mass, setting the scale using each of the string tension, Mp, or 
the S-P splitting. Here we artificially assume that the (experimentally unobserved) r]b mass 
lies 40 MeV below T assuming linear fall in the hyperfine splitting from J/ip with -p. This 
will introduces at most an 0(1/1000) error in the bare quark mass after fixing the mass of 
the spin-averaged S-state to the T system. 



1. Hyperfine Splitting 



A linear potential model argument gives the result that the quantity (M^ — M^^) is 
constant, and we compare the simulation data using M2 to define the {My + Mps) factor in 
Figure 0. Here we have included the light and heavy light data points obtained by UKQCD 



27| on the same lattices, and we use Mi for the light data points. 

Here Mp has been used to set the scale so that the data are normalised to experiment 
using Mp — M^, and the overall shape is significant. The experimental points are approx- 
imately constant for the light and heavy-light systems, while lying slightly higher for J/tp. 
This probably results from the quarkonium S-states being significantly determined by the 
Coulomb part of the potential. The simulation data are very flat above the D meson and 
the quarkonium points do not demonstrate the "emergent" behaviour shown by experiment, 
which we attribute to the Coulomb pole increasing the wavefunction at the origin. Since the 
data is normalised to experiment at the light quark end of the plot, it appears that in the 
quenched approximation the hyperfine splittings fall disproportionately quickly with quark 
mass when compared with the real world in the n, K, D sector. 

Figure presents the extrapolation of the quarkonium hyperfine splitting to charm and 
beauty using the string tension to scale both axes in the extrapolation. The mass dependence 
of the splitting is a straight line through the origin, in agreement with the non-relativistic 
expansion of QCD, and in contrast to the extrapolations obtained with the tree-level clover 
and unimproved Wilson actions, p8|,|^ , as seen in Figure The size of the hyperfine splitting 



is in clear disagreement with that in the J/ system. 



2. Fine Structure 

We compute the P-states using the 7-based operators where possible; they were found 
to be significantly statistically cleaner, as can be seen by comparing figures |^ and ^ Fur- 
thermore data existed for the degenerate quarkonium with the 7 operators for each of the 
K values. The 2"*"+ state, however, required extra derivative source propagators for one of 
the quarks in the meson; these were generated only for /t = 0.12600, so that the fine split- 
tings involving the 2+"*" state were extrapolated linearly using non-degenerate mesons, to the 
point where the non-degenerate pseudoscalar inverse kinetic mass was equal to the r]c. mass. 
We generated correlation functions for operators in both the T2 and E representations and 
found them to be degenerate within our statistical error. We take the fitted mass from the 
T2 representation which we found to be statistically slightly cleaner. 

This analysis is not ideal but errors introduced may be estimated using the non- 
degenerate Hamiltonian in [|29|. This is easily re- written in terms of the the combinations 
(mi -|- 7712) and (mi — m2) giving that the corrections to the fine structure Hamiltonian 
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for non-degenerate quarks are proportional to . In extrapolating to charm from the 

nearby k = 0.12600, the correction to the fine splitting is of order 10%, which is lower than 
the level of statistical error. Extrapolations to the beauty systems via this method however 
are problematic. 

The extrapolations are performed using degenerate quarks and the 7-based operators for 
splittings not involving the 2++ state. The mq dependence of the fine structure is shown 
in Figures [1^, [TT|, |T^ and |T3|. As can be seen, we resolve completely the x triplet, however. 



the calculated splitting between the he and Xci is consistent with zero. We obtain estimates 
for the ratio m(3Pi)-m(3Po) 1-6(6) for charmonium and the systematics in extrapolating to 
bottomonium leave the extrapolation not worth considering here. The experimental value 
is 0.48 for charmonium and 0.66 for bottomonium. The same splitting when calculated in 
NRQCD has been shown to suffer from significant discretisation effects [jl^, so that without 
calculation of the same quantity at other lattice spacings we cannot distinguish between 
discretisation and quenching effects. 



3. Spectrum Predictions 

The results obtained using each of the string tension, the p mass, and the quarkonium 
S-P splitting, extrapolated to both the J/tp and T systems can be found in Table [VIII . 



As an illustration of the resulting data we present diagrams for the spectrum obtained 



using the string tension for both the J/^jj^ Figure 0, and the T, Figure |T5|, systems. The 
spectrum obtained for J/ip is broadly comparable to experiment, with the spin splittings 
underestimated, while we expect T the results show severe quenching effects in both the 
hyperfine and S — P splittings. 

Our values for the hyperfine splitting in charmonium however, agree with the FNAL 



quenched data which has been included in reviews by Shigemitsu ||3^ and Davies |T^. For 
both systems the discrepancies in the fine structure may well include strong discretisation 
effects. It is common for the lattice spacing to be set independently in each of the J/ijj 
and T systems via the S — P splitting at that quark mass. This would certainly improve 
the spectrum diagram we obtain for T, however it is felt that, since this defines a different 
quenched theory for each system, this obscures the fact that the one quenched theory cannot 
simultaneously reproduce the spectra of these systems. Previous calculations of charmonium 
P-states with a relativistic lattice action have been made P, |3ll , ^|32| , however, this is the 
first relativistic lattice calculation that competely resolves the x triplet. 



III. SCALING BEHAVIOUR OF HYPERFINE SPLITTING 

We also performed a less extensive calculation using the tadpole- improved clover action 
at /3 = 6.2, where we consider only the S-states. The simulation used 150 sample gauge 
configurations from a quenched distribution at /? = 6.2 on a 24^ x 48 lattice. Five heavy 
quark masses corresponding to, 

ft^heavy G {0.12000,0.12300,0.12600,0.12900,0.13200}, 
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were used, all having the tadpole improved clover action with Csw = 1.44239, corresponding 
to Mo = 0.88506 taken from the plaquette. The data were fitted using a similar procedure 
to that described previously. 

With only two calculations, an extrapolation to the continuum using a linear model is 
ill-advised, and so we plot the lattice spacing dependence without performing a fit. The 
scaling behaviour of the J/ip hyperfine splitting is shown in Figure |T6|, and can be seen 
to be near scaling with both the string tension and Mp scales, in clear disagreement with 
experiment. The charmonium hyperfine results are consistent with those of El-Khadra et al 
|T[]. Based on this it is reasonable to assume that our results on the T hyperfine splitting 
underestimate the (as yet undetermined) experimental value. 

Future calculations should produce a result for the hyperfine splitting which scales using 
the quarkonium S-P splitting to set the lattice spacing. This will probably not agree with 



those using the low momentum observables since NRQCD calculations produce a result |T7 
for the S-P splitting in units of Mp which scales and is 30% lower than experiment for T; 
this induces an inconsistency between the quarkonium hyperfine splitting scaled with each 
of these quantities. The mass dependence of the quenched S-P splitting will, of course, mean 
that the level of discrepancy between the lattice spacing definitions will differ in the J/ip and 
T systems. For comparison we can estimate the value NRQCD calculations of the hyperfine 
splitting would yield using Mp to set the scale by assuming Ahyperfine oc — ; feeding the 
change in scale into both the selected quark mass and the splitting then yields a 40-50% 
reduction in the hyperfine splitting. The NRQCD quenched continuum result appears to 
be around 35 MeV Jl^ using the S-P splitting to set the scale, with a 30% uncertainty due 
to Uq ambiguities and higher order relativistic corrections. This gives an estimate for the 
central value of around 17-21 MeV if Mp were used to set the scale. 



IV. CHARM QUARK MASS 

We first fix the bare lattice mass, t^charm, by requiring the kinetic mass of the pseudoscalar 
matches that of the rjc using the string tension to set the scale. We then obtain values for 
the pole mass of the charm quark using two methods involving either mi or 1712, which 
differ due to lattice spacing effects. 

The first method evaluates the binding energy at that bare quark mass: 

aBi = aMiiQQ) - 2aM[^, (15) 
where the perturbative pole mass contains its one loop correction obtained from |]19| 



«<Q^ = E/W (16) 



1=0 



The binding energy, defined via Mi, is then subtracted from the physical meson mass to 
give the charm pole mass. 

MP°'^ = M{1S) - Bi (17) 
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The other method employed relates the one loop kinetic mass to the bare quark mass 



using the all orders in amg renormalisation constant |19 

M2 = Zm^ X m2{MD, 
where M^^ contains its one loop correction and 



m2{aMi 



aMi 



sinh aMi 
1 + sinh aMi 



(18) 



(19) 



The largest source of error in the calculation is the coupling g"^. Ideally we use the 
Brodsky-Lepage-Mackenzie scale q* for the self-energy. We have not however evaluated this 
scale and instead proceed in a similar manner to p^ , |35[] by varying the scale at which we 
evaluate av{q) over the range ^ < g < | to estimate the likely higher order error. 

The scaling behaviour obtained for the two methods for obtaining the charm pole mass 



is given in Figure The error is dominated by the uncertainty in the coupling, however it 
is plausible that a consistent continuum limit is reached by the two methods, and suggest 
the likely range for the pole mass is 1.25 GeV < m^h^'^ < 1.55 GeV. We find the pole mass 
is, naturally, rather insensitive to the upper limit of q and strongly sensitive to the lower 
limit on q due to the nature of the running coupling. 



We note that the perturbative errors are large here, in particular Gray et al find 
a large O(a^) correction in the relation between the pole and MS masses. Due to this 
uncertainty we do not convert our results to MS and instead quote only the pole mass. 
This is not to say, of course, that we believe that the pole mass is better defined, but rather 
we are reluctant to apply the one loop formula when a large two loop correction may exist. 



V. CONCLUSIONS 

The quenched hyperfine splittings in quarkonium systems scale within statistical error 
between the two calculations with either the string tension or Mp used to set the scale, 
and are in clear disagreement with experiment for J/ip and with phenomenological model 
predictions for T. Our results for the hyperfine splitting in charmonium are in agreement 
with , and those for T are in plausible agreement with NRQCD results when a correction 
for the different inverse lattice spacings used is made. Systematic errors are estimated to be 
at the 5% level. 

We make the first lattice calculation of the complete IP triplet fine structure using a 
relativistic quark action, and find that the inconsistencies of the spectrum with experiment 
in the calculation at (5 = 6.0, namely the low hyperfine splitting and the mass dependence 
of the S — P splitting are consistently explicable by an underestimation of the Coulomb 
coefficient at hadronic length scales in the quenched approximation. We give an estimate 
for the pole mass of the charm quark, with an error which is dominated by perturbative 
uncertainties. 
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FIGURES 

FIG. 1. Estimated error in mean kinetic energy estimated at tree level as a fraction of 
the total kinetic energy, A^jn at (3 = 6.0 and /3 = 6.2 due to discretisation effects. Wc used 
an approximation to potential model values {v'^) in this estimate. Estimated systematic 
errors are of the order 5%. 
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FIG. 2. Heavy-heavy pseudoscalar and vector effective masses for a meson containing 
degenerate quarks with k = 0.12600 at /3 = 6.0. These plots are considered representative 
of the rest of the data set. The hyperfine sphtting is clearly seen. The circles are points 
from local-local correlation functions, and the fancy squares come from the fuzzed source, 
local sink correlation functions. 
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FIG. 3. Quarkonium XOjXI; and h states using 7 operators for a meson containing 
degenerate quarks with k = 0.12600 at /3 = 6.0. The 0"'' — 1"*" splitting can be clearly seen. 
The 1"*"+ and 1"' states are not resolved. The circles are points from local-local correla- 
tion functions, and the fancy squares come from the fuzzed source, local sink correlation 
functions. 
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FIG. 4. Quarkonium XOjXIj X2 states using differential operators for meson 
containing degenerate quarks with k = 0.12600 at /3 = 6.0. The complete fine structure 
of the X triplet is clearly resolved, even before correlated differences are taken. 
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FIG. 5. Fit to the dispersion relation for the pseudoscalar meson containing degener- 
ate quarks with k = 0.12600 at (3 = 6.0. The fit is quadratic, and a significant difference 
between the fitted coefficient and the naive derivative obtained from the first two points 
is found, making the quadratic fit necessary. 
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FIG. 6. Plots show the S — P scaled to the string tension at /3 = 6.0, and the inverse 
lattice spacing obtained from quarkonium S — P splitting as a function of the inverse 
kinetic pseudoscalar meson mass. No significant deviations from linearity are found within 
statistical error. 
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FIG. 7. Hyperfine splitting across entire mass range at /3 = 6.0. The first three 
data points are degenerate hght quarks (open circles), there are fifteen heavy- light data 
points (fancy diamonds), and five degenerate heavy-heavy data points (fancy squares). 
The heavy-heavy points show no sign of a rise relative to heavy-light points, in contrast 
to the experimental behaviour. This shows a clear defect of the quenched approximation. 
Here we have factorised my — rripg and used a kinetic mass factor m2{PS) + m2{V) and 
a pole mass factor mi(F) — mi{PS) for all but the light quarks, for which we use solely 
mi terms. 
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FIG. 8. Extrapolation of quarkonium hyperfine splitting in the inverse pseudoscalar 
meson mass at /3 = 6.0. The axes in this plot are scaled to the string tension. The 
ratio of the hyperfine splitting to the string tension is clearly underestimated by the data 
compared with experiment. The hyperfine splitting extrapolates to zero in the static limit. 
The axes are scaled using the string tension. 
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FIG. 9. Previous UKQCD heavy quark extrapolation of quarkonium hyperfine split- 
ting using various actions. The hyperfine splitting clearly disappears before the static 
limit, in contrast to the tadpole improved data set. 
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FIG. 10. Heavy quark extrapolation of the quarkonium Xci — XcO fine splitting at 
/3 = 6.0 using degenerate quarks and the 7 operators. This splitting is clearly non-zero, 
but appears to be significantly underestimated at this lattice spacing. Whether this effect 
is due to quenching or discretisation errors cannot be understood without corresponding 
simulations on other lattice spacings. The axes are scaled using the string tension. 
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FIG. 11. Heavy quark extrapolation of quarkonium Xci — he fine splitting at /? = 6.0 
using degenerate quarks and the 7 operators. This is statistically consistent with zero at 
all masses simulated. The axes are scaled using the string tension. 
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FIG. 12. Heavy quark extrapolation of quarkonium Xc2 — he fine splitting at /3 = 6.0. 
Here we use non-degenerate quarks k\ G {0.114,0.118,0.122,0.126} and K2 = 0.126, and 
plot the difference between the mass of the 2++ state using the differential operator and the 
l"*""*" using the 7 operator versus the inverse kinetic mass of the corresponding pseudoscalar 
meson. The systematic errors in making predictions for degenerate quarkonia are expected 
to be below statistical errors for cc, but are poorly controlled for bb. The results are scaled 
using the string tension. 
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FIG. 13. Heavy quark extrapolation of quarkonium Xc2 — XcO fine splitting at /? = 6.0 
in a similar manner to figure The axes are scaled using the string tension. 




FIG. 14. Charmonium spectrum diagram using the string tension to set the scale at 
j3 = 6.0. Broad agreement is seen, but it appears that the spin structure is underesti- 
mated. 
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FIG. 15. Upsilon spectrum diagram using string tension to set scale at /3 = 6.0. 
Clearly the S—P splitting is underestimated when using the same scale as the charmonium 
simulation. This effect is universally believed to be due to quenching. It is commonly 

absorbed by using the S-P splitting to set the scale in a system dependent manner, which 
would involve a loss of predictive power here. Presenting our data in this manner shows 
that any single quenched theory cannot simultaneously get both spectra correct. 
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FIG. 16. Scaling behaviour of J/ip and T hyperfine splitting between j3 = 6.0 and 
/3 = 6.2. Results scale well with both the string tension and Mp. We have no value for 
the S — P splitting at /? = 6.2; it is possible that one would obtain results with S — P 
which scale, but yield different continuum limit (particularly for T), the discrepancy being 
a quenching error. 
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FIG. 17. Scaling Behaviour of fnF^^j.-ai- The error here is dominated by the uncertainty 
in q*. It is plausible that the results from the method 1 (diamonds) and method 2 (squares) 
will agree in the continuum limit. 
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TABLES 





TABLE I. 


6 = 6.0 heavy qiiark kappas 






aMpsiQQ) 


Fuzzing Radius 


MR Iterations 


0.13000 

0.12600 
0.12200 
0.11800 
0.11400 


0.9283(7) 
1.1618(6) 
1.3755(6) 
1.5751(6) 
1.7644(6) 


CO CO CO CO CO 


95 
75 
65 
60 
55 


TABLE II. Meson operators 


State 




Operators 


'So 

's, 

'Pi 
'Po 

-'Pi 

-'P2 


0-+ 

1— 
1+- 
0++ 

1++ 

2++ 


V'7iV' 

- 7.)Aj}V' E rep 
Vi{7i + 7j AilV' T rep 




TABLE III. 


d = 6.0 fitted timeslice ranges 


State 


Pscudoscalar 


Vector 


P-states 


Fit Range 


8-20 (FL,LL) 


8-20 (FL,LL) 


6-18 (FL,LL) 
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TABLE IV. Heavy-heavy S states 











)\/f_, 


rSeUQO \^oUuU,oUUUJ 


U.OO 


U.Do 


n Q98'^l'7^ 

u.yzoo^^ ( j 


1/11 /''^^ 


Jr^SeuaO ^^ZDUU,ZDUUJ 


U.OO 


U.Di 


i.iDiy(^o j 






n QQ 




l.O I OO^O J 




Pseudo (1800,1800) 


1.2 


0.25 


1.5751(6) 


1.97(1) 


Pseudo (1400,1400) 


1.4 


0.093 


1.7644(6) 


2.14(1) 


Vector (3000,3000) 


1.0 


0.44 


0.9780(9) 


1.49(2) 


Vector (2600,2600) 


0.95 


0.52 


1.1990(8) 


1.66(2) 


Vector (2200,2200) 


1.0 


0.44 


1.4052(7) 


1.83(1) 


Vector (1800,1800) 


1.2 


0.26 


1.5998(7) 


2.00(1) 


Vector (1400,1400) 


1.4 


0.12 


1.7854(7) 


2.17(1) 



TABLE V. Heavy-heavy P states 



State 


xVdof 


Q 


Mgs 


Mex 


0++ (3000,3000) 


1.1 


0.37 


1.161(7) 


1.76(10) 


0++ (2600,2600) 


0.92 


0.56 


1.376(7) 


1.89(7) 


0++ (2200,2200) 


0.91 


0.57 


1.578(6) 


2.04(5) 


0++ (1800,1800) 


1.0 


0.45 


1.770(6) 


2.19(4) 


0++ (1400,1400) 


1.2 


0.28 


1.955(7) 


2.34(4) 


1++ (3000,3000) 


0.62 


0.90 


1.193(7) 


1.76(7) 


1++ (2600,2600) 


0.56 


0.94 


1.404(7) 


1.91(6) 


1+" (2200.2200) 


0.65 


0.87 


1.602(6) 


2.07(5) 


1++ (1800,1800) 


0.79 


0.73 


1.791(7) 


2.22(4) 


1++ (1400,1400) 


0.91 


0.58 


1.972(7) 


2.37(4) 


1+- (3000,3000) 


0.56 


0.94 


1.193(8) 


1.67(6) 


1+- (2600,2600) 


0.42 


0.99 


1.402(7) 


1.86(5) 


1+- (2200,2200) 


0.51 


0.96 


1.602(7) 


2.04(5) 


1+- (1800,1800) 


0.70 


0.83 


1.792(7) 


2.21(4) 


1+" (1100.1100) 


0.86 


O.Gl 


1.975(7) 


2.37(4) 
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TABLE VI. Momenta or kinetic mass analyses 



Momentum Directions 





(0,0,0) 




(1,0,0) 


(0,1,0) 


(0,0,1) 


(1,1,0) 


(1,0,1) 


(0,1,1) 




(1,1,1) 




(2,0,0) 


(0,2,0) 


(0,0,2) 


TABLE VIL Heavy-heavy M2{PS) 


Stale 


x'/dof 




Pscudo (3000,3000) 


1.0 


1.024(11) 


Pseudo (2600,2600) 


0.14 


1.350(11) 


Pseudo (2200,2200) 


0.14 


1.686(12) 


Pseudo (1800,1800) 


0.6 


2.05(1) 


Pseudo (1400,1400) 


0.9 


2.43(1) 


Pseudo (2600,1400) 


0.12 


1.87(1) 


Pscudo (2600,1800) 


0.06 


1.69(1) 


Pseudo (2600,2200) 


0.025 


1.51(1) 



TABLE VIII. (3 = 6.0 heavy-heavy mass splittings 



Splitting (MeV) 




Mp 


QQ S-P Splitting 


Experiment 


J/lp - 7]c 


68(2) 


65(2) 


79(4) 


117 


^Pi-S 


418(13) 


408(16) 




458 


Xc2 - Xcl 


81(28) 


78(27) 


93(33) 


46 


Xcl - XcO 


51(7) 


49(7) 


59(11) 


95 


\c2 — XcO 


133(28) 


128(28) 


153(35) 


141 


-rih 


22(1) 


21(1) 


31(3) 




^Pi-S 


366(17) 


358(20) 




460 


Xb2 — Xbl 


43(30) 


42(30) 


56(32) 


21 


Xbl - XbO 


26(9) 


25(9) 


33(10) 


32 


Xb2 - XbO 


65(30) 


63(27) 


85(32) 


53 
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TABLE IX. /? = 6.2 heavy-heavy mass sphttings 



Splitting (MeV) 




Mp 


Experiment 




72(3) 


66(3) 


117 




22(3) 


20(3) 
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